/**
 * @file test_basic_distrs.cpp
 *
 * Unit testing of basic distributions
 * 
 * @author Dahua Lin
 */


#include <gtest/gtest.h>
#include <slipp/prob_base/bernoulli_distr.h>
#include <slipp/prob_base/uniform_distr.h>
#include <slipp/prob_base/duniform_distr.h>
#include <slipp/prob_base/normal_distr.h>
#include <slipp/prob_base/discrete_distr.h>
#include <slipp/prob_base/exponential_distr.h>

#include "eigen_test.h"
#include "random_test.h"

#include <limits>

using namespace sli;
using namespace sli::test;
typedef sli::std_rstream<sli::default_std_rand_eng> rstream_t;


/************************************************
 *
 *   Bernoulli distribution
 *
 ************************************************/

TEST( BernoulliDistr, Basics )
{
	bernoulli_distr D1;

	EXPECT_EQ( D1.p(), 0.5 );
	EXPECT_EQ( D1.q(), 0.5 );
	EXPECT_DOUBLE_EQ( D1.mean(), 0.5 );
	EXPECT_DOUBLE_EQ( D1.var(), 0.25 );
	EXPECT_DOUBLE_EQ( D1.entropy(), math::log(2.0) );

	bernoulli_distr D2(0.4);

	EXPECT_EQ( D2.p(), 0.4 );
	EXPECT_EQ( D2.q(), 0.6 );
	EXPECT_DOUBLE_EQ( D2.mean(), 0.4 );
	EXPECT_DOUBLE_EQ( D2.var(), 0.24 );
	EXPECT_DOUBLE_EQ( D2.entropy(), -( 0.4 * math::log(0.4) + 0.6 * math::log(0.6) ) );

	bernoulli_distr::param_type param(0.4);
	bernoulli_distr D3(param);

	EXPECT_EQ( D3.p(), 0.4 );
	EXPECT_EQ( D3.q(), 0.6 );
	EXPECT_TRUE( D3.param() == D2.param() );
	EXPECT_FALSE( D3.param() != D2.param() );
}


TEST( BernoulliDistr, Evaluation )
{
	bernoulli_distr D(0.4);

	typename pdf_evaluator_of<bernoulli_distr>::type ev = D.get_pdf_evaluator();

	EXPECT_DOUBLE_EQ( ev.pdf(true),  0.4 );
	EXPECT_DOUBLE_EQ( ev.pdf(false), 0.6 );

	EXPECT_DOUBLE_EQ( ev.logpdf(true),  math::log(0.4) );
	EXPECT_DOUBLE_EQ( ev.logpdf(false), math::log(0.6) );

	barow_t x(2);
	x(0) = false;
	x(1) = true;

	darow_t r_pdf = ev.pdf(x);
	darow_t r_logpdf = ev.logpdf(x);

	EXPECT_DOUBLE_EQ( r_pdf(1),  0.4 );
	EXPECT_DOUBLE_EQ( r_pdf(0), 0.6 );

	EXPECT_DOUBLE_EQ( r_logpdf(1),  math::log(0.4) );
	EXPECT_DOUBLE_EQ( r_logpdf(0), math::log(0.6) );
}


TEST( BernoulliDistr, Sampling )
{
	double p = 0.4;
	bernoulli_distr D(p);

	rstream_t rs;
	typename sampler_of<bernoulli_distr, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;

	barow_t x1(N);
	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_NEAR( double(x1.count()) / double(N), p, 1.0e-3);

	barow_t x2 = sampler.draw((size_t)N);

	EXPECT_NEAR( double(x2.count()) / double(N), p, 1.0e-3);
}



/************************************************
 *
 *   Uniform distribution
 *
 ************************************************/

TEST( UniformDistr, Basics )
{
	uniform_distr D1;

	EXPECT_EQ( D1.a(), 0.0 );
	EXPECT_EQ( D1.b(), 1.0 );
	EXPECT_EQ( D1.span(), 1.0);
	EXPECT_DOUBLE_EQ( D1.mean(), 0.5 );
	EXPECT_DOUBLE_EQ( D1.var(), 1.0 / 12 );
	EXPECT_DOUBLE_EQ( D1.entropy(), std::log(1.0) );

	uniform_distr D2(2.0, 5.0);

	EXPECT_EQ( D2.a(), 2.0 );
	EXPECT_EQ( D2.b(), 5.0 );
	EXPECT_EQ( D2.span(), 3.0 );
	EXPECT_DOUBLE_EQ( D2.mean(), 3.5 );
	EXPECT_DOUBLE_EQ( D2.var(), 0.75 );
	EXPECT_DOUBLE_EQ( D2.entropy(), std::log(3.0) );

	uniform_distr::param_type pa(2.0, 5.0);
	uniform_distr D3(pa);

	EXPECT_EQ( D3.a(), 2.0 );
	EXPECT_EQ( D3.b(), 5.0 );
	EXPECT_EQ( D3.span(), 3.0 );
	EXPECT_TRUE( D3.param() == D2.param() );
	EXPECT_FALSE( D3.param() != D2.param() );
}

TEST( UniformDistr, Evaluation )
{
	uniform_distr D(2.0, 6.0);

	darow_t x = darow_t::LinSpaced(10, 0.0, 9.0);
	darow_t r0_pdf = darow_t::Zero(10);
	r0_pdf.segment(2, 5).fill(0.25);
	darow_t r0_logpdf = darow_t::Constant(10, -std::numeric_limits<double>::infinity());
	r0_logpdf.segment(2, 5).fill(math::log(0.25));

	typename pdf_evaluator_of<uniform_distr>::type
	ev = D.get_pdf_evaluator();

	index_t n = x.size();
	darow_t r1_pdf(n);
	darow_t r1_logpdf(n);
	for (index_t i = 0; i < n; ++i)
	{
		r1_pdf(i) = ev.pdf(x(i));
		r1_logpdf(i) = ev.logpdf(x(i));
	}

	EXPECT_TRUE( test_equal(r0_pdf, r1_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r1_logpdf) );

	darow_t r2_pdf = ev.pdf(x);
	darow_t r2_logpdf = ev.logpdf(x);

	EXPECT_TRUE( test_equal(r0_pdf, r2_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r2_logpdf) );

}


TEST( UniformDistr, Sampling )
{
	double a = 2.0;
	double b = 6.0;
	uniform_distr D(a, b);

	rstream_t rs;
	typename sampler_of<uniform_distr, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;
	darow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	ASSERT_TRUE( (x1 >= a).all() );
	ASSERT_TRUE( (x1 <= b).all() );
	EXPECT_NEAR( x1.mean(), D.mean(), 1.0e-3 );
	EXPECT_NEAR( x1.minCoeff(), a, 1.0e-3 );
	EXPECT_NEAR( x1.maxCoeff(), b, 1.0e-3 );

	darow_t x2 = sampler.draw((size_t)N);

	ASSERT_TRUE( (x2 >= a).all() );
	ASSERT_TRUE( (x2 <= b).all() );
	EXPECT_NEAR( x2.mean(), D.mean(), 1.0e-3 );
	EXPECT_NEAR( x2.minCoeff(), a, 1.0e-3 );
	EXPECT_NEAR( x2.maxCoeff(), b, 1.0e-3 );
}


/************************************************
 *
 *   Discrete Uniform distribution
 *
 ************************************************/

TEST( DUniformDistr, Basics )
{
	duniform_distr<> D1;

	EXPECT_EQ( D1.a(), 0 );
	EXPECT_EQ( D1.b(), 1 );
	EXPECT_EQ( D1.span(), 2);
	EXPECT_DOUBLE_EQ( D1.mean(), 0.5 );
	EXPECT_DOUBLE_EQ( D1.var(), 0.25 );
	EXPECT_DOUBLE_EQ( D1.entropy(), std::log(2.0) );

	duniform_distr<> D2(2, 5);

	EXPECT_EQ( D2.a(), 2 );
	EXPECT_EQ( D2.b(), 5 );
	EXPECT_EQ( D2.span(), 4 );
	EXPECT_DOUBLE_EQ( D2.mean(), 3.5 );
	EXPECT_DOUBLE_EQ( D2.var(), 1.25 );
	EXPECT_DOUBLE_EQ( D2.entropy(), std::log(4.0) );

	duniform_distr<>::param_type pa(2, 5);
	duniform_distr<> D3(pa);

	EXPECT_EQ( D3.a(), 2 );
	EXPECT_EQ( D3.b(), 5 );
	EXPECT_EQ( D3.span(), 4 );
	EXPECT_TRUE( D3.param() == D2.param() );
	EXPECT_FALSE( D3.param() != D2.param() );
}


TEST( DUniformDistr, Evaluation )
{
	duniform_distr<> D(2, 5);

	index_t m = 10;
	iarow_t x(m);
	for (index_t i = 0; i < m; ++i) x(i) = int32_t(i);

	darow_t r0_pdf = darow_t::Zero(m);
	r0_pdf.segment(2, 4).fill(0.25);
	darow_t r0_logpdf = darow_t::Constant(m, -std::numeric_limits<double>::infinity());
	r0_logpdf.segment(2, 4).fill(math::log(0.25));

	typename pdf_evaluator_of<duniform_distr<> >::type
	ev = D.get_pdf_evaluator();

	darow_t r1_pdf(m);
	darow_t r1_logpdf(m);
	for (index_t i = 0; i < m; ++i)
	{
		r1_pdf(i) = ev.pdf(x(i));
		r1_logpdf(i) = ev.logpdf(x(i));
	}

	EXPECT_TRUE( test_equal(r0_pdf, r1_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r1_logpdf) );

	darow_t r2_pdf = ev.pdf(x);
	darow_t r2_logpdf = ev.logpdf(x);

	EXPECT_TRUE( test_equal(r0_pdf, r2_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r2_logpdf) );
}


TEST( DUniformDistr, Sampling )
{
	int32_t a = 2;
	int32_t b = 5;
	duniform_distr<> D(a, b);

	rstream_t rs;
	typename sampler_of<duniform_distr<>, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;
	iarow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_EQ( x1.minCoeff(), a );
	EXPECT_EQ( x1.maxCoeff(), b );
	EXPECT_NEAR( double((x1 == 2).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x1 == 3).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x1 == 4).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x1 == 5).count()) / double(N), 0.25, 1e-3);

	iarow_t x2 = sampler.draw((size_t)N);

	EXPECT_EQ( x2.minCoeff(), a );
	EXPECT_EQ( x2.maxCoeff(), b );
	EXPECT_NEAR( double((x2 == 2).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x2 == 3).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x2 == 4).count()) / double(N), 0.25, 1e-3);
	EXPECT_NEAR( double((x2 == 5).count()) / double(N), 0.25, 1e-3);

}



/************************************************
 *
 *   Normal distribution
 *
 ************************************************/

inline double nrmd_entropy(double sigma)
{
	return 0.5 * math::log((2 * SLI_M_PI * SLI_M_E) * math::sqr(sigma));
}


TEST( NormalDistr, Basics )
{
	normal_distr D1;

	EXPECT_EQ( D1.mu(), 0.0 );
	EXPECT_EQ( D1.sigma(), 1.0 );
	EXPECT_DOUBLE_EQ( D1.mean(), 0.0 );
	EXPECT_DOUBLE_EQ( D1.var(), 1.0 );
	EXPECT_DOUBLE_EQ( D1.entropy(), nrmd_entropy(1.0));

	normal_distr D2(3.0, 2.0);

	EXPECT_EQ( D2.mu(), 3.0 );
	EXPECT_EQ( D2.sigma(), 2.0 );
	EXPECT_DOUBLE_EQ( D2.mean(), 3.0 );
	EXPECT_DOUBLE_EQ( D2.var(), 4.0 );
	EXPECT_DOUBLE_EQ( D2.entropy(), nrmd_entropy(2.0));

	normal_distr::param_type pa(3.0, 2.0);
	normal_distr D3(pa);

	EXPECT_EQ( D3.mu(), 3.0 );
	EXPECT_EQ( D3.sigma(), 2.0 );
	EXPECT_TRUE( D3.param() == D2.param() );
	EXPECT_FALSE( D3.param() != D2.param() );
}


inline double nrmd_pdf(double mu, double sigma, double x)
{
	double sig2 = math::sqr(sigma);
	double c = math::rsqrt(SLI_M_PI_2X * sig2);
	double u = math::sqr(x - mu) / (2 * sig2);
	return c * math::exp(- u);
}


TEST( NormalDistr, Evaluation )
{
	double mu = 3.0;
	double sig = 2.0;
	normal_distr D(mu, sig);

	darow_t x = darow_t::LinSpaced(21, -5.0, 5.0);
	index_t n = x.size();

	darow_t r0_pdf(n);
	darow_t r0_logpdf(n);

	for (index_t i = 0; i < n; ++i)
	{
		r0_pdf(i) = nrmd_pdf(mu, sig, x(i));
		r0_logpdf(i) = math::log(r0_pdf(i));
	}

	typename pdf_evaluator_of<normal_distr>::type
	ev = D.get_pdf_evaluator();

	darow_t r1_pdf(n);
	darow_t r1_logpdf(n);
	for (index_t i = 0; i < n; ++i)
	{
		r1_pdf(i) = ev.pdf(x(i));
		r1_logpdf(i) = ev.logpdf(x(i));
	}

	EXPECT_TRUE( test_approx(r0_pdf, r1_pdf, 1.0e-12) );
	EXPECT_TRUE( test_approx(r0_logpdf, r1_logpdf, 1.0e-12) );

	darow_t r2_pdf = ev.pdf(x);
	darow_t r2_logpdf = ev.logpdf(x);

	EXPECT_TRUE( test_approx(r0_pdf, r2_pdf, 1.0e-12) );
	EXPECT_TRUE( test_approx(r0_logpdf, r2_logpdf, 1.0e-12) );
}


TEST( NormalDistr, Sampling )
{
	double mu = 3.0;
	double sig = 2.0;
	normal_distr D(mu, sig);

	rstream_t rs;
	typename sampler_of<normal_distr, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;
	darow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_NEAR( x1.mean(), D.mean(), 5.0e-3 );
	EXPECT_NEAR( (x1 - mu).square().mean(), D.var(), 2.0e-2 );

	darow_t x2 = sampler.draw((size_t)N);

	EXPECT_NEAR( x2.mean(), D.mean(), 5.0e-3 );
	EXPECT_NEAR( (x2 - mu).square().mean(), D.var(), 2.0e-2 );

}


/************************************************
 *
 *   Discrete distribution
 *
 ************************************************/

TEST( DiscreteDistr, Basics )
{
	const int K = 4;
	double p[K] = {0.1, 0.4, 0.3, 0.2};
	darow_t pv = make_arow(K, p);

	double mv = 0.4 * 1 + 0.3 * 2 + 0.2 * 3;
	double ent = -(
			0.1 * math::log(0.1) +
			0.4 * math::log(0.4) +
			0.3 * math::log(0.3) +
			0.2 * math::log(0.2) );

	discrete_distr<int> D(pv);

	ASSERT_EQ( D.K(), 4 );
	ASSERT_EQ( D.p(0), 0.1 );
	ASSERT_EQ( D.p(1), 0.4 );
	ASSERT_EQ( D.p(2), 0.3 );
	ASSERT_EQ( D.p(3), 0.2 );

	EXPECT_DOUBLE_EQ( D.mean(), mv );
	EXPECT_DOUBLE_EQ( D.entropy(), ent );

	discrete_distr<int>::param_type param(pv);

	discrete_distr<int> D2(param);

	ASSERT_EQ( D2.K(), 4 );
	ASSERT_EQ( D2.p(0), 0.1 );
	ASSERT_EQ( D2.p(1), 0.4 );
	ASSERT_EQ( D2.p(2), 0.3 );
	ASSERT_EQ( D2.p(3), 0.2 );
}

TEST( DiscreteDistr, Evaluation )
{
	const int K = 4;
	double p[K] = {0.1, 0.4, 0.3, 0.2};
	darow_t pv = make_arow(K, p);

	discrete_distr<int> D(pv);

	index_t m = 10;
	iarow_t x(m);
	for (index_t i = 0; i < m; ++i) x(i) = int32_t(i);

	darow_t r0_pdf = darow_t::Zero(m);
	for (int i = 0; i < K; ++i) r0_pdf(i) = p[i];

	darow_t r0_logpdf = darow_t::Constant(m, -std::numeric_limits<double>::infinity());
	for (int i = 0; i < K; ++i) r0_logpdf(i) = math::log(p[i]);

	typename pdf_evaluator_of<discrete_distr<int> >::type
	ev = D.get_pdf_evaluator();

	darow_t r1_pdf(m);
	darow_t r1_logpdf(m);
	for (index_t i = 0; i < m; ++i)
	{
		r1_pdf(i) = ev.pdf(x(i));
		r1_logpdf(i) = ev.logpdf(x(i));
	}

	EXPECT_TRUE( test_equal(r0_pdf, r1_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r1_logpdf) );

	darow_t r2_pdf = ev.pdf(x);
	darow_t r2_logpdf = ev.logpdf(x);

	EXPECT_TRUE( test_equal(r0_pdf, r2_pdf) );
	EXPECT_TRUE( test_equal(r0_logpdf, r2_logpdf) );
}

TEST( DiscreteDistr, Sampling )
{
	const int K = 4;
	double p[K] = {0.1, 0.4, 0.3, 0.2};
	darow_t pv = make_arow(K, p);

	discrete_distr<int> D(pv);

	rstream_t rs;
	typename sampler_of<discrete_distr<int>, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;
	darow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_EQ( x1.minCoeff(), 0 );
	EXPECT_EQ( x1.maxCoeff(), K-1 );
	EXPECT_NEAR( double((x1 == 0).count()) / double(N), p[0], 1e-3);
	EXPECT_NEAR( double((x1 == 1).count()) / double(N), p[1], 1e-3);
	EXPECT_NEAR( double((x1 == 2).count()) / double(N), p[2], 1e-3);
	EXPECT_NEAR( double((x1 == 3).count()) / double(N), p[3], 1e-3);

	iarow_t x2 = sampler.draw((size_t)N);

	EXPECT_EQ( x2.minCoeff(), 0 );
	EXPECT_EQ( x2.maxCoeff(), K-1 );
	EXPECT_NEAR( double((x2 == 0).count()) / double(N), p[0], 1e-3);
	EXPECT_NEAR( double((x2 == 1).count()) / double(N), p[1], 1e-3);
	EXPECT_NEAR( double((x2 == 2).count()) / double(N), p[2], 1e-3);
	EXPECT_NEAR( double((x2 == 3).count()) / double(N), p[3], 1e-3);
}


TEST( DiscreteDistr, SortedSampling )
{
	const int K = 4;
	double p[K] = {0.1, 0.4, 0.2, 0.3};
	darow_t pv = make_arow(K, p);

	darow_t sw(K);  sw << 0.4, 0.3, 0.2, 0.1;
	iarow_t si(K);  si << 1, 3, 2, 0;

	discrete_distr<int> D(pv);

	rstream_t rs;
	sorted_discrete_sampler<rstream_t, int> sampler = D.get_sorted_sampler(rs);

	ASSERT_TRUE( test_equal(sampler.sorted_weights(), sw) );
	ASSERT_TRUE( test_equal(sampler.sorted_values(), si) );

	index_t N = 1000000;
	darow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_EQ( x1.minCoeff(), 0 );
	EXPECT_EQ( x1.maxCoeff(), K-1 );
	EXPECT_NEAR( double((x1 == 0).count()) / double(N), p[0], 1e-3);
	EXPECT_NEAR( double((x1 == 1).count()) / double(N), p[1], 1e-3);
	EXPECT_NEAR( double((x1 == 2).count()) / double(N), p[2], 1e-3);
	EXPECT_NEAR( double((x1 == 3).count()) / double(N), p[3], 1e-3);

	iarow_t x2 = sampler.draw((size_t)N);

	EXPECT_EQ( x2.minCoeff(), 0 );
	EXPECT_EQ( x2.maxCoeff(), K-1 );
	EXPECT_NEAR( double((x2 == 0).count()) / double(N), p[0], 1e-3);
	EXPECT_NEAR( double((x2 == 1).count()) / double(N), p[1], 1e-3);
	EXPECT_NEAR( double((x2 == 2).count()) / double(N), p[2], 1e-3);
	EXPECT_NEAR( double((x2 == 3).count()) / double(N), p[3], 1e-3);
}


/************************************************
 *
 *   Exponential distribution
 *
 ************************************************/

inline double expd_entropy(double beta)
{
	return 1.0 + math::log(beta);
}


TEST( ExponentialDistr, Basics )
{
	exponential_distr D1;

	EXPECT_EQ( D1.beta(), 1.0 );
	EXPECT_EQ( D1.lambda(), 1.0 );
	EXPECT_DOUBLE_EQ( D1.mean(), 1.0 );
	EXPECT_DOUBLE_EQ( D1.var(), 1.0 );
	EXPECT_DOUBLE_EQ( D1.entropy(), expd_entropy(1.0));

	exponential_distr D2(2.0);

	EXPECT_EQ( D2.beta(), 2.0 );
	EXPECT_EQ( D2.lambda(), 0.5 );
	EXPECT_DOUBLE_EQ( D2.mean(), 2.0 );
	EXPECT_DOUBLE_EQ( D2.var(), 4.0 );
	EXPECT_DOUBLE_EQ( D2.entropy(), expd_entropy(2.0));

	exponential_distr::param_type pa(2.0);
	exponential_distr D3(pa);

	EXPECT_EQ( D3.beta(), 2.0 );
	EXPECT_TRUE( D3.param() == D2.param() );
	EXPECT_FALSE( D3.param() != D2.param() );
}


inline double expd_pdf(double lambda, double x)
{
	return x >= 0 ? lambda * math::exp(-lambda * x) : 0.0;
}


TEST( ExponentialDistr, Evaluation )
{
	double beta = 2.0;
	double lambda = 1.0 / beta;
	exponential_distr D(beta);

	darow_t x = darow_t::LinSpaced(13, -1.0, 5.0);
	index_t n = x.size();

	darow_t r0_pdf(n);
	darow_t r0_logpdf(n);

	for (index_t i = 0; i < n; ++i)
	{
		r0_pdf(i) = expd_pdf(lambda, x(i));
		r0_logpdf(i) = math::log(r0_pdf(i));
	}

	typename pdf_evaluator_of<exponential_distr>::type
	ev = D.get_pdf_evaluator();

	darow_t r1_pdf(n);
	darow_t r1_logpdf(n);
	for (index_t i = 0; i < n; ++i)
	{
		r1_pdf(i) = ev.pdf(x(i));
		r1_logpdf(i) = ev.logpdf(x(i));
	}

	EXPECT_TRUE( test_approx(r0_pdf, r1_pdf, 1.0e-12) );
	EXPECT_TRUE( test_approx(r0_logpdf, r1_logpdf, 1.0e-12) );

	darow_t r2_pdf = ev.pdf(x);
	darow_t r2_logpdf = ev.logpdf(x);

	EXPECT_TRUE( test_approx(r0_pdf, r2_pdf, 1.0e-12) );
	EXPECT_TRUE( test_approx(r0_logpdf, r2_logpdf, 1.0e-12) );
}


TEST( ExponentialDistr, Sampling )
{
	double beta = 2.0;
	exponential_distr D(beta);

	rstream_t rs;
	typename sampler_of<exponential_distr, rstream_t>::type sampler = D.get_sampler(rs);

	index_t N = 1000000;
	darow_t x1(N);

	for (index_t i = 0; i < N; ++i)
	{
		x1(i) = sampler.draw();
	}

	EXPECT_NEAR( x1.mean(), D.mean(), 5.0e-3 );
	EXPECT_NEAR( (x1 - beta).square().mean(), D.var(), 5.0e-2 );

	darow_t x2 = sampler.draw((size_t)N);

	EXPECT_NEAR( x2.mean(), D.mean(), 5.0e-3 );
	EXPECT_NEAR( (x2 - beta).square().mean(), D.var(), 5.0e-2 );

}





